We introduce a novel approach, which integrates unfitted finite element methods (FEMs) into finite element interpolated neural networks (FEINNs), to solve partial differential equations (PDEs) on complex geometries. In this method, the neural network (NN) is interpolated on a trial finite element (FE) space defined on the background mesh elements that intersect or lie within the physical domain. The Dirichlet boundary conditions are weakly imposed by Nitsche’s method. The method is robust to variations in Nitsche coefficients and small cut cell issues. We demonstrate its effectiveness in solving both linear and nonlinear PDEs across various 2D and 3D complex geometries, including those defined by level-set functions and Stereolithography (STL) meshes. The trained NN solutions outperform FE interpolations of the analytic solution on the same trial space, achieving several orders of magnitude improvement in H1 errors for smooth solutions, while their interpolations maintain expected h- and p-convergence rates. Furthermore, preconditioning based on the dual residual norms, applied to standard or aggregated test spaces, remarkably accelerates NN training convergence and enhances robustness to the choice of Nitsche coefficients.